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The numerical relativity session at GR18 was dominated by physics results on binary black hole 
mergers. Several groups can now simulate these from a time when the post-Newtonian equations of 
motion are still applicable, through several orbits and the merger to the ringdown phase, obtaining 
plausible gravitational waves at infinity, and showing some evidence of convergence with resolution. 
The results of different groups roughly agree. This new- won confidence has been used by these 
groups to begin mapping out the (finite-dimensional) initial data space of the problem, with a 
particular focus on the effect of black hole spins, and the acceleration by gravitational wave recoil to 
hundreds of km/ s of the final merged black hole. Other work was presented on a variety of topics, 
such as evolutions with matter, extreme mass ratio inspirals, and technical issues such as gauge 
choices. 



I. INTRODUCTION 



In the numerical simulation of comparable-mass black 
hole binary mergers, after a decade of struggle with un- 
stable codes, several groups are now obtaining reliable 
gravitational wave signals, and are competing closely to 
investigate the parameter space. 



This topic dominated the plenary talk by B. 
Briigmann, as well as about half of the talks submitted to 
the parallel session B2 Numerical Methods. To celebrate 
the fact that binary black hole simulations are at last 
obtaining astrophysics results, a session was held jointly 
with the Bl Relativistic Astrophysics parallel session. 7 
contributed talks were given in B1/B2 and 24 more in 
B2. There were 6 posters in B2. 



Much work is also being done on numerical simula- 
tions of astrophysical scenarios involving matter in gen- 
eral relativity, in particular core collapse and neutron star 
binary mergers, but at GR18 this interesting and active 
area was represented by few contributed talks. There 
were no talks on computer algebra. 



Given the remarkable sudden progress in binary black 
hole evolutions since 2005, and the remarkable similar- 
ity of the results from different groups, I begin with an 
overview of this field. I have aimed to make the astro- 
physical and historical overviews (Secs. lnAl and llTB|) ac- 
cessible to researchers outside numerical relativity, and 
a description of the state of the art (Sec. Ill C[) useful 
to beginning graduate students in the field. After that, 
the summaries of individual talks are likely to be too 
technical to be accessible to any but the specialists, and 
too short to tell the specialists anything new. There- 
fore, where I could identify a paper or e-print related to 
a specific talk, I give that reference. 



II. BLACK HOLE BINARY MERGERS 
A. Astrophysical background 

Many stars are in binary systems, and many of these 
are expected to consist, at the end of their life, of 
two black holes, two neutron stars, or a neutron star 
and a black hole. All binary systems lose orbital en- 
ergy through emitting gravitational radiation, and merge 
eventually, although the time scale on which this happens 
depends on their initial separation and masses. Compact 
objects such as neutron stars and black holes approach 
each other very closely before merger, which means that 
they can emit significant amounts of gravitational radia- 
tion just before merger, releasing an energy up to a few 
percent of the total mass of the system. 

These gravitational waves are very much harder to de- 
tect than the same energy in light, basically because the 
frequencies are very much lower. This has the twin conse- 
quences that the instrument has to be very much larger, 
and that, unlike a CCD camera for light, it cannot de- 
tect individual quanta. For these reasons, no gravita- 
tional waves have yet been measured directly. (How- 
ever, the time series of several pulsars in close binaries 
agrees very accurately with the energy loss through grav- 
itational waves predicted by general relativity.) 

Mergers of compact object binaries are expected to 
be the principal source of gravitational waves with fre- 
quencies of tens to thousands of Hertz. The large inter- 
ferometric detectors LIGO, GEO, VIRGO and TAMA 
are now observing in the hundreds of Hertz frequency; 
they were reviewed in the plenary talk of S. Whitcomb. 
LIGO is funded for an upgrade that should allow it to 
see binary mergers in the local group; this should give it 
enough sources to see something. (See also the plenary 
talks by D. Shadock on space-based detectors, and by M. 
A. Papa on gravitational wave astronomy.) 

In the mid 1990s, the binary black hole problem seemed 
the natural problem for numerical relativity to tackle as 
soon as computers were big enough to allow simulations 
in 3D: it is entirely described by smooth solutions of the 
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Einstein equations, without the need to model matter. 
While the two black holes are still far apart each is ap- 
proximately a Kerr black hole, and is parameterised by 
only its mass and spin, and their orbit is parameterised by 
only its ellipticity and size. The gravitational field itself 
has an infinite number of degrees of freedom, but after a 
few orbits the gravitational wave field is effectively deter- 
mined by the history of the orbital motion. This means 
that the parameter space of binary mergers which start 
sufficiently far apart and where accretion into the black 
holes can be neglected is effectively finite-dimensional. 



B. History of numerical simulations 

Ten years and perhaps 500 person-years later, the 
problem had not been solved because the numerical simu- 
lations were stopped by numerical instabilities that could 
not be overcome by just using more resources. In hind- 
sight, there were overlapping problems concerning the 
mathematical formulation which interacted with, and 
were sometimes confused with, problems of numerical 
discretisation. These were mainly of two types: the use of 
ill-posed formulations of the initial-boundary value prob- 
lem for general relativity, and gauge choices not appro- 
priate to the symmetries of the problem. The evolution 
time before the code broke down was gradually extended, 
with the first full orbit achieved by Brugmann in 1997 [l|, 
but only a very few orbits became possible until 2005. 

A first stable simulation with an essentially unlimited 
number of orbits up to and through the merger was an- 
nounced by F. Pretorius in March 2005 2]. He solved 
the Einstein equations in modified harmonic coordinates. 
Their principal part is then a set of wave equations with 
characteristics on the physical lightconcs. However, some 
lower-order terms were modified in order to modify the 
harmonic time slicing. Other key elements of his success 
include, but were not limited to, imposing outer bound- 
ary conditions by compactifying the numerical domain at 
spatial infinity, the excision of a spacetime region inside 
each black hole (singularity excision), the use of adap- 
tive mesh refinement, and damping the constraints by 
friction-like lower-order terms. 

In November 2005, an unlimited number of orbits plus 
merger was achieved again with very different methods, 
at the same time but independently by groups at the Uni- 
versity of Texas/Brownsville [3( and at NASA/Goddard 
0. These groups used the "BSSN" formulation of the 
Einstein equations [5] , which had already been successful 
in neutron star and core collapse simulations (see for ex- 
ample @ ) . The key additional ingredient for binary black 
hole simulations was a new version of an old method for 
dealing with the singular interior of the black holes called 
the "puncture method" Q , where inside each black hole 
is a wormhole, rather than a collapsed star. 



C. State of the art 

These two formulations of the problem are still the 
dominant ones today with only minor changes. The 
BSSN formulation is now always used together with a 
"Bona-Masso 1+log" type slicing condition and a "hy- 
perbolic T-driver" type shift condition. With this gauge, 
BSSN is hyperbolic, although some of its characteristics 
are spacelike and some timelike with respect to the phys- 
ical light cones. There seems to be a trend to gauge 
drivers with time derivatives in the direction normal to 
the slice, rather than along the lines of constant spatial 
coordinates. For example, (d t + (3 l di)a = —2aK, rather 
than dtot = —2aK, and similarly for the evolution of the 
shift. This avoids a breakdown of hyperbolicity when 
gauge speeds coincide, and allows the final black hole to 
settle down to a Killing coordinate system. 

Encouraged by Pretorius' success, other groups are 
now experimenting with modified harmonic coordinates 
both for vacuum and matter evolutions, and his work 
seems to have influenced even the evolution of the lin- 
earised Einstein equations for modelling extreme mass 
ratio inspirals, such as stellar mass black holes falling 
into supermassive black holes. 

The most reliable initial data now seem to be produced 
by converting a snapshot of a post-Newtonian (PN) orbit 
into "puncture" type initial data for full general relativity 
by solving the Hamiltonian and momentum constraint. 
This PN orbit is typically constructed by calculating the 
adiabatic shrinking of a conservative PN orbit. To reduce 
eccentricity, the post-Newtonian equations of motion can 
be evolved over hundreds of orbits until the eccentricity 
has been radiated away. It is puzzling that the conversion 
works so well given that this also involves a translation 
of gauge from the post-Newtonian to the full general rel- 
ativity framework. 

Outer boundaries remain primitive - some in common 
use are not known to lead to well-posed initial-boundary 
value problems, and are known not to be compatible with 
the constraint equations. Nevertheless, numerical bound- 
ary conditions exist that are empirically stable, and the 
strategy of moving them as far out as possible, typi- 
cally by using nested boxes of Cartesian coordinate grids, 
seems to work. The location of the outer boundary does 
not seem to seriously affect the gravitational wave sig- 
nals. 

Waves are extracted on an approximately spherical 
surface either by finding the Zerilli gauge-invariant lin- 
ear perturbation with respect to a Schwarzschild back- 
ground, or by constructing an appropriate tetrad and 
calculating the Newman-Penrose scalar ^4. Here it does 
seem to be important to push this sphere as far out as 
possible. 

Black hole masses and spins are calculated by find- 
ing a dynamical horizon, and constructing approximate 
Killing vectors for use in Komar-type integrals, by en- 
ergy and angular momentum balance arguments, or by 
fitting quasinormal ringdowns to the known ringdown of 
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Kerr. Approximate ADM quantities integrated over a 
large sphere also seem to agree with these diagnostics. 



D. Contributed talks 

The contributed talks were dominated by physics re- 
sults which appear to agree between research groups. 
Most talks focused on the effect of the spin of the two 
black holes on gravitational recoil ( "kicks" ) . It appears 
that velocity of several hundred km/ s are generic, and 
several thousand km/s are possible. This is astrophysi- 
cally relevant as it may completely eject the merged black 
hole from its galaxy. Note that as velocity is naturally 
measured in units of the speed of light, and vacuum grav- 
ity is scale-invariant, these velocities are invariant under 
an overall scaling of the initial data. Therefore the re- 
sults reported here would apply equally to stellar mass 
and supermassive black hole binaries, as long as the two 
masses are of the same order of magnitude. 

P. Diener [9j reported on joint work by the AEI and 
LSU groups. Initial data were quasi-circular, determined 
by an effective potential method. Kicks up to 115km/ s 
were achieved for non-spinning unequal mass binaries 
(mass ratio 0.36), while kicks up to AAOkm/s were ob- 
tained with equal mass black holes with spins anti-aligned 
parallel to the orbital angular momentum. The post- 
Newtonian analysis of gravitational recoil by Kidder [8| 
was found to be qualitatively valid beyond its domain of 
applicability. 

M. Hannam reviewed the work of the Jena group. They 
had concentrated on a configuration with equal masses, 
and spins anti-aligned in the orbital plane ( "superkicks" ) 
[io| . The magnitude depended sinusoidally on the angle 
of the spins in the orbital plane. Because of the high 
symmetry, most of the energy was radiated as I = 2, m = 
±2 waves, and the kick could be estimated through the 
energy difference between these. The velocities obtained 
were of the order of 2500km/ s (for J = 0.723A/ 2 ), in 
agreement with what had been found by the Brownsville 
group [|| (which was not represented at GR18). 

Work was also under way on a bank of numerical tem- 
plate wave forms for non-spinning mergers . Full nu- 
merical evolutions based on post-Newtonian initial data 
could be compared to a continued post-Newtonian evo- 
lution for up to 9 more orbits before the merger, and 
agreed very accurately. Going to sixth order accuracy in 
the spatial finite differencing had been crucial to obtain 
a small enough phase error (2 degrees over 9 orbits) for 
this. Moreover, very accurate "hybrid" wave forms could 
be constructed by matching a post-Newtonian inspiral 
to black hole ringing, using information from numerical 
simulations to match them. This might be used to fill 
the parameter space with templates. 

P. Laguna reported work from the PSU group 
on spins parallel to the orbital angular momentum (up 
to AOOkm/s) and superkick configurations, and stressed 
that the Kidder post-Newtonian formula for spin-orbit 



interactions, although out of the domain of its validity, 
could be used as a heuristic formula with a small num- 
ber of free parameters to be determined by full numeri- 
cal simulations. Quasi-circular orbits were compared to 
post-Newtonian ones over 9 orbits, and agreed, although 
to lower precision than that obtained by the Jena group. 
The moving punctures method was robust. Waveform 
comparisons between groups were now needed. 

B. Kelly presented the work of the Goddard group. 
Hybrid waveforms agreed accurately with fully numeri- 
cal ones. The code was fifth-order accurate. The effects 
of unequal masses and spins were being analysed. He 
described an analysis of the gravitational wave recoil in 
multipoles [131 ] . with the force seen as a sum of mainly 
three products of pairs of multipoles, and stressed that 
the integral of momentum over time was not monotonic 
( "antikicks" ) . This model could also explain why spins in 
the orbital plane produced larger kicks. However, accre- 
tion would tend to align spins with orbital angular mo- 
mentum. Finally, the effective one-body approach was 
a step forward in covering the parameter space with ap- 
proximate wave forms. 

M. Scheel [T3| presented work by the Caltech group 
that uses a modified harmonic formulation. (In con- 
trast to Pretorius, their formulation is reduced to first 
order in space and time, and consistent boundary con- 
ditions are applied at finite radius. In contrast to all 
other groups except the Meudon group they use spectral 
methods rather than finite differencing in space.) Cur- 
rently their code experienced coordinate problems during 
merger, but spectral methods allowed for the currently 
most accurate inspiral simulations in full general relativ- 
ity. Post-Newtonian based initial data (equal mass, no 
spin, zero eccentricity) evolved for 30 orbits until merger 
agree accurately for the first 15 orbits, and then signifi- 
cant disagreement could be shown, after taking into ac- 
count the different gauges of what is being compared. 
It was found that quasi-circular initial data are actually 
slightly eccentric. 

F. Pretorius [H| was interested in unusual threshold 
behaviour in binary black hole merger, rather than as- 
trophysical wave forms. For non-spinning black holes of 
comparable mass, at the threshold of immediate merger, 
prompt merger was delayed and replaced by a whirling 
phase in which the number of orbits n depends on the 
impact parameter b as expn ~ \b — 6*| 7 . Although fully 
nonlinear, this behaviour was also shown by point par- 
ticle geodesies in Kerr. During the whirl, 1-1.5% of the 
energy was radiated per orbit. He also offered an expla- 
nation of superkicks in terms of each black hole expe- 
riencing frame dragging by the other. He expected no 
more surprises in generic orbits, except perhaps in the 
limit of extreme rotation, or infinite boost. 

D. Pollney ^ presented work of the AEI and LSU 
group on spin interactions. Distinctive features in the 
wave forms could be associated with the spins. The 
momentum flux could be estimated in particular by the 
terms Qt 2 Qt-3 and Q^Qsii- 
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D. Shoemaker [T7j presented a quantitative study by 
the PSU group of the influence of spurious radiation in 
the initial data on the merger, by artificially adding I = 
m = 2 Teukolsky waves inside the binary. This changed 
the ADM mass but not the angular momentum, with the 
merger and ringdown relatively unaffected. 

I. Hinder from the PSU group focused on the com- 
parison with post-Newtonian evolutions over 9 orbits to 
merger. Quasi-circular puncture initial data actually 
showed slight eccentricity. Disagreement with PN evo- 
lutions became more noticeable in the last 3 orbits. The 
gravitational waves matched the post-Newtonian ones 
well, although the eccentricity showed, as well as a de- 
pendence on resolution towards the end. The wave form 
still converged to 3rd to 4th order until the merger. 

S. Husa [181 ] from the Jena group spoke on a variety of 
topics arising in evolutions. Higher order finite differenc- 
ing was needed for accuracy and, surprisingly, seemed to 
work well with punctures, which are not smooth. 4th or- 
der Kreiss-Oliger dissipation was important, as well as 
using an advection stencil for the shift terms. Phase 
error converged to 6th order, but not amplitude error; 
this could be fixed by re-parameterising by phase. Very 
low eccentricity in the initial data (obtained by evolv- 
ing the post-Newtonian equations over many orbits) was 
needed to see any disagreement between full numerical 
and 3.5PN evolutions. 

J. Gonzalez, also at Jena, reported on unequal mass 
mergers, with a mass ratio of up to 10. Gauge problems 
arose which may be related to the damping term in the 
shift driver. The kick speeds agreed well with [li|. As the 
mass ratio increased, less energy was radiated in I = 2, 
and more in / = 3, 4, 5. 

W. Tichy [2(| reported on work at FAU on binaries 
with J = 0.8A/ 2 , resulting in kicks of up to 2500km/ s. 
With 10 levels of mesh refinement, evolutions with an 
outer boundary at 240M could be run on a 32Gb work- 
station. The largest error seemed to come from the ex- 
traction radius. 



III. CONTRIBUTED TALKS ON OTHER 
TOPICS 

A. Physics results 

R. de Pietri [2l[ described simulations of instabilities 
in rapidly rotating stars with a T-law equation of state 
(initially polytropic, and no shocks form). The thresh- 
old of instability was well approximated by a Newtonian 
analysis. A bar (m = 2 dominated) appeared initially, 
but later resolved into an m = 1 or m = 3 dominated 
structure. 

S. Liebling described the application of a parallel adap- 
tive mesh refinement code in harmonic coordinates to 
different physics problems. Binary boson stars [22] could 
usefully be compared with binary black holes to see how 
the gravitational waves emitted depend on the internal 



structure. Here, the relative phase of the two stars pro- 
vides an additional arbitrary parameter. Neutron star 
binaries had also been simulated [23]. The old problem 
of the critical collapse of Brill waves had been taken up 
again, although the critical amplitude had only been ap- 
proached to within 1%. Subcritical scaling gave 7 ~ 0.23, 
A ~ 0.75, which is quite different from the original re- 
sults [24| 7 — 0.36, A ~ 0.60. Magneto-hydrodynamics 
had been implemented [25[ . 

U. Sperhake [26| reported simulations of head-on col- 
lisions of black holes which had different internal struc- 
ture, namely Brill-Lindquist (BL) initial data (each black 
hole is a wormhole to a separate internal infinity), Mis- 
ner data (both black holes connect to the same internal 
infinity) and approximate initial data obtained by super- 
posing two Kerr-Schild (KS) slices through Schwarzschild 
(where the slices end at the future spacelike singularity) . 
(His is the only BSSN code that can stably move ex- 
cised black holes, which is necessary for evolving these, 
non-puncture initial data. Pretorius also uses moving ex- 
cision, but in the harmonic formulation). KS had much 
more spurious initial radiation. BS and Misner agreed 
closely, but there was a small discrepancy between these 
and KS in the wave form which might be due to numer- 
ical error. 

B. Kol reviewed a number of problems in higher dimen- 
sions, often string-inspired, which might usefully be in- 
vestigated by numerical relativity In particular he spec- 
ulated on a possible relation between Choptuik's critical 
solution in the gravitational collapse of a scalar field in 
4D, and the pinching off of a black string in Br~ x,x x S 1 

m 



B. Gauge and boundaries 

The talks by D. Brown [H and D. Garfinkle [H| were 
closely related to the talk of S. Husa [3(| given in parallel 
session A3: Why does the "moving punctures" approach 
to representing black holes work in numerical relativity? 
The initial data 0] represent a wormhole opening out 
to another universe in the Kruskal solution. However, 
evolving with the slicing condition now used by most bi- 
nary black hole groups ("Bona-Masso 1+log"), numerical 
error allowed the solution inside the black hole to jump 
to a slicing of the Kruskal spacetime where the wormhole 
is asymptotically cylindrical and ends at the internal fu- 
ture null infinity, rather than at the internal spacelike 
infinity. The slicing could then become Killing both in- 
side and outside the black hole. At the same time a 'T- 
driver" type shift condition removed grid points from the 
interior, producing in effect excision by under-resolution. 

There were two other talks on coordinate choices for 
black hole evolutions. D. Hilditch [3lJ discussed slicings 
of asymptotically flat spacetimcs that are asymptotically 
null combined with radial compactification, such that the 
speed of outgoing waves remains constant. Numerical 
tests of the spherical wave equation on Schwarzschild re- 
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solved outgoing waves to arbitrarily large radii on a small 
number of grid points. L. Lindblom discussed generali- 
sations of the "harmonic coordinate driver" coordinate 
conditions used by Pretorius. 

J. Seiler presented an implementation by the AEI 
group of boundary conditions for the generalised har- 
monic formulation of the Einstein equations which are 
at once compatible with the constraint and result in a 
well-posed initial-boundary value problem (32j. These 
had been implemented using summation by parts tech- 
niques for finite differencing and the code was now being 
tested. 

L. Brewin [33| reviewed his work on smooth lattice nu- 
merical relativity: in contrast to Regge calculus, space- 
time is not locally flat, and the legs of the simplices are 
geodesies. 

C. Evolution codes 

J. Thornburg review ongoing work on an adaptive 
mesh refinement code in double null coordinates for the 
linearised vacuum Einstein equations in harmonic gauge. 
The code has been tested with the scalar self-force on cir- 
cular orbits in Schwarzschild [34[. It is to be applied to 
the calculation of the self-force on generic particle orbits 
in the Schwarzschild, and later Kerr spacetime. 

B. Zink reviewed progress on the application of high 
order, multipatch, summation-by-parts finite differenc- 
ing techniques [3a ] to the simulation of spacetimes with 
perfect fluid and magnetohydrodynamics matter. Initial 
applications are to a rotating black hole with an accretion 
disk. A particular strength of the multipatch approach 
was that it allowed the combination of high resolution 
in the radial direction with low resolution in the angular 
directions. 

J. Novak [3f| reviewed a fully constrained nonlinear 
3D evolution code where only two variables are evolved: 
in a suitable linearised limit these represent the two po- 
larisations of gravitational waves. The spatial gauge is 
Dirac gauge, defined with respect to a flat background 
metric. The code uses spherical coordinates and spec- 
tral methods. Boundary conditions are imposed at finite 
distance which are absorbing for quadrupolar waves on 
Minkowski. It was planned to evolve a single black hole 
with horizon boundary conditions, and to include matter. 

D. Initial data 

J. Read [37j | reported on numerical work in the he- 
lical Killing field approximation to the binary inspiral 



problem, which could be seen as an alternative to the 
waveless (conformally flat) approximation. A scalar field 
toy model and a neutron star binary had been imple- 
mented. The equations are mixed hyperbolic-elliptic. 
Convergence at present required the outer boundary to 
be very close in, but this might be improved by an outer 
patch where dthij = is imposed instead. 

B. Kelly [38| presented initial data of the puncture 
type in full general relativity but based on a snapshot of 
a post-Newtonian binary orbit, and with approximately 
correct gravitational wave content also based on that or- 
bit. The constraints are currently not solved, but are 
violated only at 0[(v/c) 5 ]. 



E. Analysis 

J. L. Jaramillo described the application of a re- 
fined Penrose inequality in axisymmetry, namely A < 
8tt(M 2 + V-A/ 4 + J 2 ), as a geometric test for a new ap- 
parent horizon finder (40| . and for the outermost char- 
acter of marginal trapped surfaces. Numerical evidence 
of Dain's proposal for characterizing Kerr data by the 
saturation of such inequality was also presented. 

G. Cook [4l} presented methods for finding an approx- 
imate rotational Killing vector in numerical evolutions, 
which could then be used in a Komar-type integral to 
compute the spin of a black hole. These are constructed 
by minimising an appropriated measure of the residue 
in the Killing equation, resulting in a system of coupled 
elliptic equations. This was being tested, and might pro- 
vide a new tool for analysing the horizon geometry. 

N. Bishop presented progress on the extraction of grav- 
itational waves at scri by matching to an outgoing null 
grid. Problems with numerical noise had been over- 
come, and the extraction had been tested in the head- 
on collision of two black holes in a harmonic code with 
constraint-preserving boundary conditions. Full Cauchy- 
charactcristic matching now seemed feasible. 

A. Norton [42j presented animations of a spatial co- 
ordinate system that rotates with a rotating star in a 
central region but is fixed at infinity, without getting in- 
creasingly entangled: it is based on the double cover of 
50(3) by S 3 , and so returns to its original state after 
two rotations. 
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